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Abstract 

We propose a numerical method for computing all eigenvalues (and 
the corresponding eigenvectors) of a nonlinear holomorphic eigenvalue 
problem that lie within a given contour in the complex plane. The 
method uses complex integrals of the resolvent operator, applied to 
at least k column vectors, where k is the number of eigenvalues inside 
the contour. The theorem of Keldysh is employed to show that the 
original nonlinear eigenvalue problem reduces to a linear eigenvalue 
problem of dimension k. No initial approximations of eigenvalues and 
eigenvectors are needed. The method is particularly suitable for mod- 
erately large eigenvalue problems where k is much smaller than the 
matrix dimension. We also give an extension of the method to the 
case where k is larger than the matrix dimension. The quadrature 
errors caused by the trapezoid sum are discussed for the case of an- 
alytic closed contours. Using well known techniques it is shown that 
the error decays exponentially with an exponent given by the product 
of the number of quadrature points and the minimal distance of the 
eigenvalues to the contour. 

1 Introduction 

We consider nonlinear eigenvalue problems of the form 

T{z)v = 0, V eC"',v y^O,z eQ, {1] 
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where T : i7 — )■ C"*''" is assumed to be holomorphic in some domain C C. 
The computation of all eigenvalues and eigenvectors inside fl usually requires 
to solve two problems (see [I3],[T] for recent reviews) : 

1. Approximate localization and separation of eigenvalues in suitable do- 
mains resp. intervals, 

2. accurate computation of eigenvalues and associated eigenvectors by an 
iterative method. 

The global problem of localization can be substantially simplified if minimum- 
maximum characterizations similar to the linear case hold [23] , [20] . Voss and 
co-workers have combined these principles with locally convergent methods of 
Arnoldi or Jacobi-Davidson type (see [21], [2], [22]), and in this way provided 
an effective means for computing all eigenvalues. 

Another case where both problems can be solved, is for polynomials 

p 

T{z) = J2 - zoY, Tj e C™'™. 

j=0 

This eigenvalue problem can be reduced to a linear eigenvalue problem of 
dimension pm, and this is the path taken by the MATLAB routine polyeig. 
Quite a few papers in the literature either analyze this linearization approach 
or generalize methods from linear eigenvalues to the polynomial case. 

In the general holomorphic case we just have a power series near each 

oo 

T{z) = J2 - ZoY, \z-Zo\ small, Tj E C'"''". 

j=0 

One may then use polynomial truncation and the polynomial solver for get- 
ting good initial estimates of the eigenvalues. However, the success of this 
method strongly depends on the radius of convergence and on the decay of 
the coefficient matrices. Also, it may be necessary to compute power series 
at many different points in Q. 

Finally, we refer to the recent approach of Kressner [TT], who uses the 
fact that any holomorphic matrix function can be written as 

with holomorphic functions /j : i— )■ C (such a representation always exists 
for some p < rn?). Then a Newton-type iteration is devised in pjj that allows 
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to compute a group of eigenvalues and an associated subspace. Though the 
convergence of this method is surprisingly robust to the choice of initial 
values, it remains a method for solving the local problem. 

In this paper we tackle the global problem by using contour integrals, 
which seem to be the only available tool in the general holomorphic case. The 
idea is to use the theorem of Keldysh [S],[in], which provides an expansion 
of T{z)^^ in a neighborhood U G Q of a.n eigenvalue A G as follows: 



T{z)-'=J2Sj{z-\y, zeU\{\}, S-, gC'"'"^, ^-.^0. (2) 



More specifically, Keldysh' theorem gives a representation of the singular part 
in ([2]) in terms of generalized eigenvectors of T{z) and its adjoint T^{z). A 
good reference for the underlying theory is [15] which we briefly review in 
Section [2l 

Numerical methods based on contour integrals seem not to have attracted 
much attention in the past. A notable exception are exponential integrators 
and, more recently, approaches to compute analytic functions of matrices via 
suitably transformed contour integrals, see |7j,^ 13.3.2]. In particular,the 
exponential convergence of the trapezoid sum is proved in [7]. 

Our goal is to compute all eigenvalues and the associated eigenvectors that 
lie within a given closed contour F in Q. The main algorithm is described 
in Section [31 Suppose that k < m eigenvalues of ([1]) lie inside T. Then our 
method reduces the nonlinear eigenvalue problem to a linear one of dimension 
k by evaluating the contour integrals 



Here V G C™'*'^ is generally taken as a random matrix. The contour integrals 
in ([3]) are calculated approximately by the trapezoid sum. If N quadrature 
points are used, this requires to solve Nk linear systems, which is the main 
numerical effort. As a consequence, our method is limited to moderately 
large nonlinear eigenvalue problems for which a fast (sparse) direct solver is 
available. 

In Section S] we apply the algorithm to several examples, showing that a 
moderate number of quadrature nodes (A^ ~ 25) is usually sufficient to get 
good estimates of eigenvalues and eigenvectors. Based on [4J, we prove in 
Section H] that the quadrature error decays exponentially with an exponent 
that depends on the product of the number of quadrature nodes and the 
smallest distance of the eigenvalues to the contour. 



oo 



K 




(3) 
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In the final Section [S] we deal with two problems that are typical for non- 
linear eigenvalue problems and that do not occur in the linear case: First, 
there can be much more eigenvalues than the matrix dimension (e.g. char- 
acteristic functions for delay equations) and, second, eigenvectors belonging 
to different eigenvalues can be linearly dependent, even if the number of 
eigenvalues is less than the matrix dimension. In Section [S] we extend our 
integral method such that it applies to the case k > m and that it can also 
handle rank defects of eigenspaces. For the extended integral method it is 
necessary to evaluate Ap from ([3]) for indices < p < 2[^] — 1. Numerical 
examples show that this extension is suitable for solving both aforementioned 
problems. 

Acknowledgement: The author thanks Ingwar Petersen for the support 
with the numerical experiments. 

2 Nonlinear eigenvalues and Keldysh' Theo- 
rem 

The material in this section is largely based on the monograph [15]. It 
contains a general study of meromorphic operator functions that have values 
in spaces of Fredholm operators of index 0. For our purposes it is sufficient 
to consider matrix valued mappings 



that are holomorphic in some open domain Q. We write this as T G H{Q, C™'"*) 
For a matrix A we denote by R{A) and N[A) its range and nuUspace, re- 
spectively. 

Definition 2.1. A number A G is called an eigenvalue of T(-) if T{X)v = 
for some v G C", v ^ 0. The vector v is then called a (right) eigenvector. By 
(t(T) we denote the set of all eigenvalues and by p(T) = Q\ cr(T) we denote 
the resolvent set. 

The eigenvalue A is called simple if 



Theorem 2.2. Every eigenvalue A G a(T) of T G if(l],C™''") is isolated, 
i.e. hi \ {A} C p(T) for some neigborhood U of A. 

Moreover, T{z) is meromorphic at A, i.e. there exist k G N and Sj G C"*'"* 
for j > —K such that 7^ and 



A(T(A)) = span{t;}, v^O T'{\)v ^ i?(T(A)). 



oo 




(4) 
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Remark 2.3. The number k is uniquely determined and called the order of 
the pole at A. 

The Theorem of Keldysh (see Theorem 12.61 below) gives a representation of 
the singular part 

-1 

in terms of (generalized) eigenvectors of T and . It goes back to Keldysh [H] 
with a proof given in . Generalizations of Keldysh' theorem were derived 
by Trofimov [19], who introduced the concept of root polynomials, and by 
Marcus and Sigal [I2] and Gohberg and Sigal [5J who used factorizations 
of operator functions. A simple direct proof was found by Mennicken and 
MoUer [T3j who later gave a concise approach to the whole theory in |15j . 

For the motivation of the algorithm in the next section it is instructive 
to first state Keldysh' theorem for simple eigenvalues. In this case Definition 
12.11 implies for the adjoint T^{z) 

N{T^{\)) = span{w} for some ^x; G C™, w ^ 0, 

w"T\\)v ^ 0. 

Without loss of generality we can normalize v and w such that 

w"T'{X)v = 1. (5) 
Then we are still free to further normalize either \w\ = 1 or |f| = 1. 

Theorem 2.4. Assume A G ^7 is a simple eigenvalue of T G H{Q, C"*'™) 
with eigenvectors normalized as in Then there is a neighborhood U dVt 
of A and a holomorphic function R G H{U, C"^'™) such that 

T{z)-^ = ^—vw^ + R{z), zeU\{\]. (6) 
z — A 

Moreover, let C C ^7 be a compact subset that contains only simple eigenval- 
ues A„, n = 1, . . . , A; with eigenvectors w„ satisfying 

T(A>„ = 0, <r(A„) = 0, <T'(A>„ = 1. (7) 

Then there is a neighborhood W of C in and a holomorphic function R G 
HiJA, C™'™) such that 

k 

^(^)"' = E -^^nw^ + ^(^)' ^ e W \ {Ai, . . . , A,}. (8) 
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Proof. The first part is a special case of Theorem 12.61 below. For the second 
part, note that eigenvalues are isolated and hence we can choose a neighbor- 
hood C <ZU <Z VL such that a{T) nU = {Ai, . . . , A^}. Then the function 

k 

is holomorphic in W fl p(T) and by the first part it is also holomorphic in 
suitable neighborhoods of A„, n = 1, . . . , fc. □ 

Definition 2.5. Let T G H{Q, C""'"") and \ E Q. 

(i) A function v G H{Q, C") is called a root function ofT at X if 

v{X) 0, T{X)v{X) = 0. 

The order of the zero z = X of T{z)v{z) is called the multiplicity of v 
at X and denoted by s{v). 

(ii) A tuple {vq, . . . ,f„-i) G (C")",?! > 1 is called a chain of generalized 
eigenvectors ( CGE) ofT at X if v{z) = X]j=o is a root function 
of T at A of multiplicity s{v) > n. 

(iii) For a given vq G N{T{X)), vq ^ the number 

r{vo) = max{s(f ) : t> is a root function of T at A with v{X) = Vq} 
is finite and called the rank of vq- 

(iv) A system of vectors in C™ 

V = {v^j,0 < j < mi - 1,1 < i < L) 

is called a canonical system of generalized eigenvectors (CSGE) ofT at 
X if the following conditions hold: 

(a) The vectors fg, . . . , form a basis of N{T{X)), 

(b) The tuple {vq, . . . , f^^.i) is a CGE of T at A for £ = 1, . . . , L, 

(c) m£ = max{r(t>o) : Vq G A^(T(A)) \ spanjiiQ : < u < £}} 
for £= 

One can show that a CSGE always exists and that the numbers are 
ordered according to 

nil > > • • • > ^L- 

They are called the partial multiplicities of T at X. With these notions we 
can state the following general theorem, see fl5\ Theorem 1.6.5]. 
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Theorem 2.6 (Keldysh). Let T e i/(fi,C'"'™) be given with p(T) 7^ 0. For 
A G a(T) let 

V = {v^j,0 < j < me - 1,1 < £ < L) 
be a CSGE of T at A. Then there exists a CSGE 

W = {w^j, 0<j<mi-l,l<£<L) 

of at A, a neighborhood W of A and a function _R G -f/^(W, C™^''") such that 

L mi mi ~j 

^W = EE(^-^)'' + ^eW\{A}. (9) 

£=1 j=l i/=0 

The system W, for which ([9]) holds, is the unique CSGE of at A that 
satisfies the following conditions 

r(wo) = rnc 

j mu 

J2 Yl = SuiSoj, 0<j <me-l,l<i,u<L, (10) 

a=0 13=1 

where 

T,- = 1t(^)(A), j>0. (11) 

Remark 2.7. Rather than using generalized eigenvectors one can also write 
T(z)^^ in terms of left and right root functions, see flbl Th.1.5.4]. 

The representation iQ shows that the order k of the pole in (j3]) is given 

by 

K = max{m^ : i = 1, . . . , L}. 

Further, the number L = dim(A^(T(A))) is the geometric multiplicity while 
Yle=i is algebraic multiplicity of A. In the semi-simple case mi = 
1,1 = 1, . . . , L, equations and (ITU]) simplify to 

L 

£=1 

wi''T\\)vl = 5,,, l<i,iy<L, 
which in case L = 1 further simplify to (j6]) and ([5]). 

Consider now all eigenvalues inside a compact set C C ^2. In the same way 
as ([8]) followed from ([6]), we obtain from Theorem I2.6l the following corollary. 
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Corollary 2.8. Let C C Q be compact and T e C'^-'"). Then C 

contains at most finitely many eigenvalues A„,n = l,...,n{C) with corre- 
sponding CSGEs 

Vn = (vf'', < j < nii^n - 1, 1 < ^ < ^n) , u = 1, . . . , n{C). 

Let 

Wn= (u^-'",0<j<m,,„-l,l<£<L„), n = l,...,n{C) 
be the corresponding CSGEs of such that 

and with = j^T'^^\Xn) 

j rni,^„ 

5Z 5Z ^j-a^'^+f^^n V'X,„-P = ^ul^Oj. < J < mi^n - 1, 1 < Z/ < 
a=0 /3=1 

Then there exists a neighborhood C dU d Vt and a function i? G H{U^ ^m,m-j 
such that for all 2; G W \ {Ai, . . . , Xn{c)} 

n{C) L„ i^l,n mi_„-j 
n=l ^=1 jr=l i/=0 

As a consequence of the corollary it follows that the order of the pole in 
(jl]) is given by 

K = max{m£,„ '■ < i < Ln, I ^ n < n{C)}. 

Consider now a contour F C i.e. a simple closed curve that has its 
interior int(r) in Q. An easy consequence of the residue theorem is the 
following result. 

Theorem 2.9. Let T G H{Q, C™'™') have no eigenvalues on the contour 
F C and denote by Xn,n = l,...,n{r) the eigenvalues in the interior 
int(F) C Q. Then with the CSGEs from Corollary 12.81 we have for any 

feH{nx) 

1- / /(.)T(.)-d. = E E E 7-n¥ E (12) 

n=l 1=1 j=l u=0 
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If all eigenvalues are simple the formula reads 

1 r "^^^ 

— / f{z)Tiz)-'dz = J2 fi^n)VnW^, (13) 

n=l 

where Vn, Wn are left and right eigenvectors corresponding to A^, and normal- 
ized according to 

w7;fT'(A„K = 1, n = l,...,n(r). (14) 

Proof. Corollary 12.81 applies to C = int(r)ur, where the function f{z)T{z)~'^ 
has residues at Xj given by the right-hand side of ffT^ . The special case 
Ln = 1, mon = 1, n = 1, . . . , n(r) yields equation (fT3|) . □ 



3 The algorithm for a few eigenvalues 

In the following we set up an algorithm for computing all eigenvalues of 
T G H{Q, C™'''") inside a given contour F in Q. We assume that the sum of 
all algebraic multiplicities 

n(r) L„ 

fc = ^^m£,n (15) 

n=l 1=1 

is less than or equal to the system dimension m. For the opposite case we 
refer to Section [51 In high-dimensional problems we actually expect to have 
k <^m. 



3.1 Simple eigenvalues inside the contour 



As in the second part of Theorem 12.91 let us assume that all eigenvalues 
Ai, . . . , A.„(r) in int(F) are simple so that k = n{r). We introduce the matrices 

V= {vi...v,),W={wi...Wk) eC^^K 

We assume that we have chosen a matrix 

V e C™^ k<l<m, 

such that 

W"V G C'^'' has rank k. (16) 

In particular, this implies rank(iy) = k. In the applications we choose V at 
random (see Section H]), so that f lT6|) can be expected to hold in a generic sense 
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if rank(iy) = k. We note that (in contrast to linear eigenvalue problems) it is 
easy to construct nonlinear eigenvalue problems for which W is rank deficient. 
However, this seems to be a nongeneric situation for typical applications. In 
addition to (|T6|) we assume 

rank{V) = k, (17) 

which again is expected to hold in generic cases. 
Next we compute the two integrals 




The evaluation of these integrals by quadrature rules is by far the most 
expensive part of the algorithm and will be discussed below. Note also, that 
in the linear case T{z) = zl — A the matrix Aq is obtained by applying to V 
the Riesz projector onto the invariant subspace associated with all eigenvalues 
inside T. 

By ([13]) we obtain 

k 

Ao = J2 ^nW^V = VW^V. (20) 

n=l 

Similarly, 

k 

Ai = Y,^nVnwlV = VAW''V, A = diag(A„,r2 = l,...,fc). (21) 

n=l 

In the next step we compute the singular value decomposition (SVD) of 
Aq in reduced form 

VW^V = Ao = VoJ^oW^ (22) 

where Vo € C'"'^ Eq = diag(ai, . . . , a,), l^o e C''^ ^V^, = 4, W^^^Wo = h- 
Note that the rank conditions (fT6l) . f|T7|) show that rank(Ao) = /c, hence Aq 
has singular values 

(Ti>...(Tfc>0 = (Tfc+i = ... = Cr;. 

By the rank condition (fT7|) we have 

R{Aq) = R{V) = R{Vo). 
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Since both, Vq and V are m x k matrices and Vq has orthonormal columns, 
we obtain 

V = VqS, S = V^V G C'^''^ nonsingular. (23) 
With ([20]), (123]) we find VqSW^V = K)SoWo^ and thus 

This relation is used to eliminate W^V from Ai = VqSAW^V . We obtain 

Vo"A^ = SAW^V = SAS-^^oW^, 

which upon multiplication by 1^0^(7^ from the right finally gives 

SAS-^ = Vo"AiWoJ:q\ (24) 

Note that the right-hand side is a computable matrix which is diagonalizable 
and has as eigenvalues exactly the eigenvalues of T inside the contour. We 
summarize the result in a theorem. 

Theorem 3.1. Suppose that T G H{Q,C^'^) has only simple eigenvalues 
Ai, . . . , Afc inside the contour F in 1] with left and right eigenvectors normal- 
ized as in f|T^ . Moreover, let a matrix V G C"'' be given such that k < I < m 
and the rank conditions f|T6]) . f|T7]) are satisfied. Then the matrix 

B = VfAiWo^Q^ G C'=•^ (25) 

given by f lTS]) .f lT^ and the SVD fl22]) . is diagonalizable with eigenvalues 
Ai, . . . , Afc. From the eigenvectors Si, . . . , G C'^ of one obtains the eigen- 
vectors of T through 

Vn = VoSn, n= l,...,k. 
Remarks 3.2. (a) For reasons of numerical stability we may replace Ai by 

ii = ^ J^iz - zo)T{z)-'Vdz = - zoAo. 

For example, in case of a circle F, one can take zo as its center. Then 
(12T]) holds with A — instead of A and the matrix B = Vf^ AiWq'Eq^ has 
eigenvalues A„ — ^o- Therefore, the eigenvalues of T are found by adding zq 
to the eigenvalues of B. 

(b) The rank conditons in the theorem are crucial. Assume, for example, 
that Aq = VW^V has rank ko < k. Then the SVD ([22]) holds with matrices 
Wo G C^'''",Vo G C'"''^'' and Eq = diag(ai, . . . , a^J. Moreover, we have 
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S e C^«'^ in 



B e C'^o-'^o in 



Finally we find B = SAS where 



S = W^VWqTiq^ satisfies SS = Ik^- Except for the case, when 5* has some 
zero columns this does not lead to a useful relation between the eigenvalues 
of B and A. For numerical computations we therefore recommend to test the 
residuals ||T(An)un||, see Section [3l3l A general cure of this rank deficient 
case is provided by the generalized algorithm in Section O which, however, is 
computationally more expensive. 



3.2 Multiple eigenvalues inside the contour 

Let us consider the general case where T G H{Q, C"^'™') has no eigenvalues on 
the contour F but may have multiple eigenvalues inside. We apply Corollary 
12.81 to the compact set C = F Uint(F) and assume that the matrix composed 
of all CSGEs that belong to eigenvalues inside F, 



V 



(t^ < J < m,,„ - 1, 1 < £ < L„, 1 < n < n(F) 



(26) 



has rank k, cf. f|T5l) . Then, using Theorem 12.91 with f{z) = 1 shows that Aq, 
as defined in flTSl) . satisfies 



£ n t,nH -fr 



n=l e=l u=0 

Further, we assume that the matrix 
has maximum rank k, where 



(27) 



W 



wi" _i_^,0<z/<m£,„-l,l<£<L„,l<n<n(F)) G C'"•^ 



•mi, 



is normalized as in Theorem 12. 6[ With Theorem 12.91 we then find 



(2^ 



n(r) L„ 



71=1 £=1 



"l«,„-l 



u=0 



u=0 



vj w^^ _i_„+ > V,: w: 



V = VAW^V, 



where A has Jordan normal form 







\ 




^ Jn,l 








1 \ 


A = 






1 Jfi 














\ 


Jn{V)) 




\ 


Jn,Ln ) 




\ 





(29) 
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As in Section [XT] the next steps are the SVD fl22]) for Aq and the computation 
of B = Vf^ AiWqTiq^ G C'^'^. Then B has eigenvalues Ai, . . . , A„(r) and its 
Jordan normal form has the same partial multiplicities as T{z). 

Theorem 3.3. Suppose that T G H{Q, C"''") has no eigenvalues on the 
contour F in and pairwise distinct eigenvalues A„,n = 1, . . . ,n(T) inside 
r with partial multiplicities rrii^n > • • • > f^L„,m n = 1, . . . , n{T). Moreover, 
assume that the matrix of generalized eigenvectors from (l26i) and the matrix 
W^V from (1271) have rank k with k given by (fT5|) . Then the matrix B G C'^''^ 
from (1251) has Jordan normal form (129]) with the same eigenvalues A„ and 
partial multiplicities „ = 1, . . . , L„, n = 1, . . . , Ti(r)). Suitable CSGEs 
for T can be obtained from corresponding CSGEs " for B via 

t;-'" = Vosf\ < J < m,,„ - 1, 1 < £ < L„, 1 < n < niV). 

Remark 3.4. Essentially, the theorem reduces the nonlinear problem for 
eigenvalues inside a contour to a linear eigenvalue problem for a x fc-matrix. 
The linear eigenvalue problem inherits the multiplicity structure of the non- 
linear problem. As usual, computing the Jordan normal form is not a stable 
process and other forms, such as the Schur form, are recommended. A closer 
look at the derivation of the algorithm (l22i) . (l24l) shows that it is sufficient to 
have a rank revealing Q-R-decomposition. One would then replace Wq^o ^ 
(I24p by the inverse of the maximum rank upper triangular submatrix. 



3.3 Quadrature and numerical realization 

The major step in the algorithm consists in evaluating the integrals (ITS]) and 
([T^ by numerical quadrature and by solving the linear systems involved in 
the evaluation of the integrand. We assume that V has a 27r-periodic smooth 
parameterization 

G C^(M,C), i^{t + 2-k) = ^{t) VtGM. 

Of particular interest is the real analytic case ip G C"^(M, C). Taking equidis- 
tant nodes tj = j = 0, . . . , N and using the trapezoid sum, we find the 
following approximations 

N-i (30) 
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where we used (p{to) = (p{tN). Similarly, 

Ai ^ A,,^ = ^ E nv{t,))-'Vv{t,)v'{t,). (31) 

In order to compute Aq^^ we need to solve Nl linear systems with different 
matrices T{ip(tj)),j = 0, . . . , A^ — 1 and with / different right-hand sides each. 
Note that we can use the solutions of these hnear systems to compute Ai^at 
at almost no extra cost. For the special case of a circle (p{t) = fi + i?e** we 
obtain the formulas 

Ao,N=§J^T{^{t,)rVexp{^-^), 

j=0 

Ai,N =^iAo,N + ^ E T(v^(t,))-Vexp(^). 
i=o 

The algorithm can be summarized as follows: 
Integral algorithm 1 

Step 1: Choose an index / < m and a matrix V G C™'' at random. 

Step 2: Compute Aq^nAi,n from (150]), CT. 

Step 3: Compute the SVD Aq^n = VEW^, where 

V e C'"'^ W e C'-', V^V = W^W = IuT. = diag((7i, as, • • • , (yi)- 

Step 4: Perform a rank test for S, i.e. find < A; < Z such that 

CTl > . . . > fXfc > tolrank > CTfe+l ^ • • • ^ CT, ^ 0. 

If A; = / then increase I and go to Step 1. 

Else let Vo = 1^(1 : m, 1 : k), Wq = W{1 ■.l,l:k) and 

Eo = diag(cri,a2, . . . ,(Tfc). 

Step 5: Compute B = Vq" Ai^nWoJ:^^ E C^'^ 

Step 6: Solve the eigenvalue problem for B 
BS = SA, S = (si . . . Sfc), A = diag(Ai, . . . , A^). 

If ||T(Aj)fj|| < tolres and A^ G int(r) accept Vj = VoSj as eigenvector and Xj 
as eigenvalue. 

Remarks 3.5. (a) If we find k = I positive singular values in Step 4 then we 
take this as an indication that there may be more than / eigenvalues (includ- 
ing multiplicities) inside F. We then increase / until a rank drop is detected 
in Step 4. 

(b) In general, it is more efficient to compute Ai^m in Step 5, when the index 
k has been determined. Then one has to store the solutions of the linear 
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systems solved during the evaluation of Aq^^. 

(c) As noted in Remark 13.2( b) the algorithm may fail due to linear de- 
pendency of (generalized) eigenvectors. Therefore, we include a test of the 
residual. Moreover, as the experiments in Section H] show, eigenvalues close 
to the contour, either inside or outside F, may lead to difficulties in the rank 
test. Therefore, the trivial test Xj G int(r) is included in Step 6 as well. 

(d) In Step 6 we assumed that eigenvalues are simple. If multiplicities occur 
or B is only brought into upper triangular form, then the eigenvalues can 
still be read off from the diagonal, and the structure of eigenvectors can be 
retrieved from VoS. 

4 Error analysis and numerical examples 
4.1 Error analysis 

Standard results on the trapezoid sum for holomorphic periodic integrands 
imply exponential convergence at a rate that depends on the number of nodes 
times the width of the horizontal strip of holomorphy, see [3],|H 4.6.5]. Ap- 
plications of these results to the computation of matrix functions via contour 
integrals appear in [7]. 

Theorem 4.1. Let / G H{S{d^,d^),C) be 27r-periodic on the strip 

S{d_,d^) = {z eC: -d_ < Imz < d± > 0. 

Then the error of the trapezoid sum 

•^0 i=0 

satisfies for all < r_ < (i_ , < r+ < (i+ 

\E^{f)\< max 1/(^)1 G(e-^''+)+ max |/(;.)| G(e-^^-), 

Im(2)=r-|, Im(2)=r_ 

where G{x) = x ^ 1. 

Remark 4.2. Note that Theorem 14.11 is a slight variation of [U 4.6.5] since 
/ is not assumed to be real on [0, 27r] and the strip S{d-, d+) can be unsym- 
metric, in general. 

In the following we state and prove the corresponding result for integrals 
over circles which will be used in the sequel. 
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Theorem 4.3. Let / G if (y4(a_, a_|_), C) be holomorphic on the annulus 

1 \z\ 

A{a_, a+) = {z : — <~F7< > 1) 

for some R > 0. Then the error of the trapezoid sum 

^^^^^ = 2^/ ^^ = exp(^), (32) 

satisfies for all 1 < p_ < a_, 1 < < a+ 

\E^{f)\ < max |/(z)| G(p;^))+ max |/(^)| G(p:^). (33) 

\z\=p+H p-\z\=H 

Proof. We use the Laurent expansion of / (see e.g. P) 

oo „ 



(34) 



which converges uniformly on compact subdomains of the annulus. By a 
simple computation, 

E (,'') = S k + l = iN,ieZ\{0}, 

^ otherwise. 

Applying to flM|l leads to 

oo 

^iv(/) = - $^(/£7vi?'^ + /-.ivi?-'^). (35) 
e=i 

From Cauchy's Theorem and a standard estimate we obtain 



R 



< ^27c p+Rmax\,\=p^R\f{z)\{p+Ry 



max 



In a similar way, 

|/-.^i?-'^|< max 1/(^)1 p:^^. 
Using these estimates in (135!) completes the proof. □ 
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The proof shows that the p_-term can be discarded in fl55]) if the principal 
term in the Laurent expansion vanishes (i.e. /fc = for k < —1). Likewise, 
the p+-term disappears when fk = for A; > 0. For the function 

f{z) = {z-X)-\ J>1, (36) 

the principal term vanishes for |A| > i? while the secondary term vanishes 
for |A| < R. Example fl36|) is crucial for the application to the meromorphic 
functions from Section [31 Therefore, we note the following explicit formula. 

Lemma 4.4. The error of the trapezoid sum (l32i) for the function (136|) in 
case > J is given as follows. 

In particular, 

Er,{{z-xrn = < > : (3^ 




Remark 4.5. If / G H{A{a_,a^),C) is meromorphic on an open neigh- 
borhood of the closed annulus A{a-,a+Y, then the estimate (133|) can be 
sharpened as follows 

E^(/) = o(a;^ + a:^). 

In order to see this, first consider the singular part that belongs to poles on 
the boundary of y4(a_, a+), and use Lemma [4.41 Then apply Theorem 14.31 to 
the remaining part on a slightly larger annulus. 

Consider a general contour F in with 27r-periodic parametrization if{t),t G 
[0, 27r]. Moreover, assume that has a 27r-periodic holomorphic extension to 
a strip 

if e H{S{d_,d+),n), ip{z + 27i) = ip{z). (39) 
For definiteness, we also assume that 

/ e int(F), < lm{z) < d+, 
' \ i int(F), < Im(2) < 0. ^ ^ 

Common examples are circles v^(2;) = -Zo + -Re'^ with z G C and ellipses 
>~p{^z) = acos(z) + bsm{z) with | Im(z)| < artanh(min(|, ^)). 
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Let g G HiVL, C), then the error of the trapezoid sum for f{z) = g{{p{z))(p'{z), 
z G 5'((i„, is 

From Theorem 14.11 we obtain an estimate 

\EM\ < $(r+)G(e-^'-+) + <l>(r_)G(e-^-), (42) 

where < r_ < (i„,0 < r+ < (i+ and $(r) = maxim(2)=.,. (7(^9(2;)) |. 

The following lemma gives a rough estimate of the right-hand sides for the 
pole function g{z) = {z — \)~^ , X E Q. 

Lemma 4.6. Let Q be bounded and let if satisfy conditions ( 1391) . fHUl) . 
Then there exist constants Ci,C2,C3 > (depending on ip, j but not on 
or A G i7) such that for dist(A, F) < C3, 

|^iv((- - A)"^)| < Cidist(A,F)-^ exp(-C2Ardist(A,F)). (43) 

Proof. For a fixed < g < 1 there are bounds |v5'(-2)| < M+ for < lm{z) < 
qd+ and \ip'{z)\ < M_ for < - Im(z) < qd_. Let C3 = max(M+rf+, M_rf_) 
and define r_|_ = Then there exists some Zj^ = s+ + zr+, < s+ < 27r 

such that 

min \X-<^{z)\= |A - (^(z+)| > |A - (^(s+)| - |(^(s+) - (^(2;+)| 

Im(2;)=r+ 

> dist(A, F) - M+r+ = (1 - g)dist(A, F). 

The first term in fH2|) can be estimated as follows 

|$(r+)|G(e-^'^+) < M+maxi^,=,^ \{y,{z) - A)-^-|G(e-^^+) 

< C(l - g)-W+dist(A, F)-^' exp (^-A^dist(A, F)-^^) . 

The second term is treated analogously. □ 

As a consequence of Lemmas 14.41 and 14.61 we obtain an exponential esti- 
mate for the errors in (130|) and (13T|) . 

Theorem 4.7. Let T G i/(i7,C) have maximum order k of poles for the 
inverse in Q, cf. Theorem 12.21 Further, let F be a simple closed contour in 
Q with cr{T) fl F = and such that the parametrization (p satisfies ( 139|) and 
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(j^0|) . Then there exist constants Ci, C2 > (depending on T and V but not 
on A^) such that the matrices from fl5Ul) . flMl) satisfy 

\\Ap - Ap,^\\ < Cid(T)-'^e-^^^^(^), p = 0, 1, 

where d{T) = mmx^a(T) dist(A, F) and d{T) = 1 if cr(T) = 0. If F is a circle 
with parametrization (p{t) = zq + i?e**, then the following estimate holds 

\\Ap - A,M\ < Ci [p^-'^+i + p^^+'^-i] , p = 0,l, 

where 

|A — zol R 
p„ = max — , p. = max r. 

Aea(T),|A-2o|<H R \ea{T),\\-zo\>R \X — Zq\ 

Combining these estimates with the well-known perturbation theory for 
singular value decompositions [IS] we find that the integral algorithm detects 
the correct rank k of Aq^^ if is sufficiently large. Further, the perturbation 
theory for simple eigenvalues [18] leads to the following corollary. 

Corollary 4.8. Let the assumptions of Theorem 13. II and of Theorem 14. 71 be 
satisfied. Let Ai, . . . , be the eigenvalues of T inside F and let Ai^at, . . . , X^^n 
be the eigenvalues from step 6 of the integral algorithm. With the notation 
from Theorem 14.71 we then have the error estimates 

max I A,- - A,,^| < Cirf(r)-"e-^^^''(^\ 

j=l,...,n{r) 

in case of a general curve satisfying f l39|) . fHOl) . and 

max \Xj - Xj,n\ < C [p^-''+^ + p^'+"-'l 

j=l,...,n{T) 

in case of a circle with radius R and center zq. 
4.2 Numerical examples 

Example 4.9. For the first test we choose a real quadratic polynomial 

T{z) = To + zTi + z^Ts, e M^°'^°, j = 0, 1, 2, (44) 

where Tq,Ti,T2 are taken at random {rand from MATLAB). In this case we 
can compare with the spectrum Upoiyeig resulting from MATLAB 's polyeig. 

Figure [Deleft) shows the result from polyeig (open circles) and the eigen- 
values from Integral algorithm 1 (filled boxes) for the data 

ip{t) = Re'\ t e [0, 2tx], R = 0.33, toWk = 10"^ toU = 10"^ (45) 
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The eight eigenvalues inside the circle are detected and well approximated 
by the integral algorithm. Figure [T] (right) shows the errors 

e{Xj) = min{|Aj - fi\ : fi e o-poiydg} 

for two characteristic eigenvalues inside the circle. Both show exponential 
decay with respect to at approximately the same rate. 
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Figure 1: Example 14. 91 Eigenvalues of a quadratic eigenvalue problem from 
polyeig (open circles) and Integral algorithm 1 (filled squares) with = 150 
(left). Difference e(Aj) of eigenvalues Ai ~ 0.30578 (filled circles) and A2 ~ 
0.0961 — 0.1315i (open circles) between polyeig and the integral algorithm 
versus the number of nodes N (right). 



|Q» . °°ooo° 



8080 

„Ho80q 



Ooo88o„ 

.o°o 8 



°°°° 



Oooo8@ 

O ° O o o 
O O O n 



0000880 



20 40 60 80 100 120 140 1' 



'N 



Figure 2: Example 14.91 Singular values versus for a fixed number of / = 11 
columns in the integral algorithm (left), reduction of the number of singular 
values by the rank test of the adaptive algorithm versus (right). 

While Figure [1] (left) results from the integral algorithm with an adaptive 
number / of columns (which yields / = 8 at A^ = 150), the computations in 
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Figure (Upright) are done with a fixed number of Z = 11 columns. For tfiis 
case we sliow the behavior of the 11 largest singular values of Aq^^ in Figure 
[2] (left). Sufficient separation of singular values already occurs at values 
N ^ 25, much smaller than 150. Figure [2] (right) shows how the adaptive 
algorithm reduces the number of singular values from / = 23 at iV = 20 to 
/ = 8 for > 95. 



Example 4.10. For the next experiment we take random complex entries 
in (jH]), a fixed number / = 10 of columns, and the same circle as in (H5ll . 
Again, the 6 eigenvalues inside the circle from polyeig are well approximated 
by the integral algorithm, see Figure [3] (left). 
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Figure 3: Example 14.101 Eigenvalues from polyeig (open circles) and eigen- 
values from the integral algorithm for a random quadratic complex matrix 
polynomial (left), singular values of integral algorithm with / = 10 columns 
versus the number of quadrature nodes for the same example (right). 

But this time the singular values do not separate as well as in Figure |2] 
(left). Two of them decay rather slowly, while two others, due to eigenval- 
ues very close but outside the contour, remain of order one. However, this 
behavior does not result in spurious eigenvalues. On the contrary, if we keep 
/ = 10 for the eigenvalue computation, then this yields the 6 eigenvalues 
inside and in addition the four eigenvalues lying closest to the contour, but 
outside. Such a behavior is also suggested by our error analyis in Section 
14.11 according to which the principle error term depends on the distance of 
eigenvalues to the contour, both for eigenvalues inside and outside. Compu- 
tational experience shows that only very small singular values (~ 10~^°) lead 
to spurious eigenvalues and these can be easily avoided by the residual test 
in Step 6. 
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Example 4.11. This example, taken from [17j and [TT], is a finite element 
discretization of a nonlinear boundary eigenvalue problem 



-u"{x) = Xu{x), < a; < 1, u{0) = = u'{l) + 



X 



A- 1 



-u(l] 



The matrix function is T{z) = Ti + jz^emC^ — zT^, where 



Ti 



m 



/2 



2 -1 
-1 1/ 



1 

6m 



/4 1 
1 ■•• 

V 



4 1 
1 2/ 



We use m = 400 and compute five eigenvalues in the interval [2,298]. 
Again Figure H] (left) shows the real eigenvalues in the circle which agree with 
those from [TT]. Note that we avoided the singularity of T at 2; = 1. The 
residuals of the computed eigenvectors and eigenvalues decay exponentially 
as expected, see Figure HI but not as smooth as in the previous examples. 
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Figure 4: Example 14. 11[ Eigenvalues from the integral algorithm for the finite 
element discretization of a nonlinear boundary eigenvalue problem (left), 
decay of residuals res(Aj) = ||T(Aj)(fj)|| for Ai ~ 24 (open circles), A2 ~ 123 
(filled circles) versus the number N of quadrature nodes for the same example 
(right). 
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Example 4.12. Consider the quadratic polynomial 



T{z) =To + {z-a){b~ z)Ti, a<beM, To,T^eR^^'^\ (46) 

where Tq has zeroes in the first column. All other entries of Tq, Ti are chosen 
at random. Then T{z) has different eigenvalues a and b with the same 
eigenvector G M"*. This is a critical case since the rank condition (fTTj) is 
violated. In Figure [5] (left) we show the results of polyeig and of the integral 
algorithm (with / = 5 and the data from (H5i) ). There are three eigenvalues 
inside the circle. Both eigenvalues a = —0.2 and b = 0.1 are missed by 
the integral method, while the third one is found, though at lower accuracy 
than in the previous examples. Figure [5] shows that only one singular value 
stays of order one when is increased. This example will be reconsidered 
in Section O 
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Figure 5: Example 14.111 Eigenvalues from polyeig (open circles) and eigen- 
values from the integral algorithm for a quadratic matrix polynomial with 
rank defect (left), singular values of integral algorithm with / = 5 columns 
versus the number of quadrature nodes for the same example (right). 



5 The algorithm for many eigenvalues 

In this section we show how the method from Section [3] can be extended to 
nonlinear eigenvalue problems with more eigenvalues than the dimension of 
the system, i.e. m < k, and to the rank deficient cases, see Remark 13.21 and 
Example KT2\ 
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5.1 Construction of algorithm 

In case m < k condition (|T7|) is always violated and there is no matrix 

V satisfying ( IT6|) . Therefore, we compute more integrals of type ( TTSj) . (|T9|) . 
namely 

Ap = — [ z^TizY^Vdz e C™'', peN. 

Here we assume that V G C"'' with I < m. In fact, in case k > m we set 

V = Im instead of making a random choice. 

From Theorem 12.91 we obtain 

Ap = VAPW^V, PEN, (47) 
and (1251) and A has the normal form 



where V,W E C"^'^ are given by 

Now we choose K E N, K > 1 and form the Km x Kl matrices 



/ ^0 



Ak-i\ 



\Ak-i ■ ■ ■ ) 
From (147|) we find the representations 

/ V \ 



(A, 
\Ak 



Ak \ 
A2K-1 J 



(48) 



5n 



and 



/ V \ 



(49) 



Bi 



A{W"V ■■■ A^-^W"V) 



(50) 



yvA^'-'j 



We assume that K has been chosen such that the following rank condition 
holds 

/ ^ \ 



rank 



k. 



(51) 



\VA^'^J 

The smallest index having this property is called the minimality index in 
In case k > m this can be expected to hold if we choose 



{K - l)m < k < Km. 



24 



In case k < m with rank(l^) < k (see Remark 13.2( b)) the following lemma 
shows that fl^T]) holds for K larger than the sum of the maximal ranks at all 
eigenvalues. 

Lemma 5.1. Let the assumptions of Corollary 12.81 be satisfied. Then the 
rank conditon flSTjl holds with k as defined in f|T5l) for 

n(C) 

K >} max rrii 

^ 1<£<L„ 
n=l 



Proof. Let M„ = maxi<£<£,„ mi^n and M = Yln=i Assume that VK^ 
0,j = 0,...,M — 1 for some x G C™. For any n G {!,..., n(C)} and 
< /3 < M„ — 1 consider the polynomial 

n{C) 

By our assumption = V^P„,/3(A)x. We partition according to (l29l) 

F= (^1 •■■ K(c)) , K = ■■■ K,L„) , K,, = (^;^'" ■■■ t;^';^^!) 

Using {Jfi — \n)^^^ = for n 7^ n we obtain 

L„ n(C) 

= Vn,i n (^-'^ - ^")*'' ("^-'^ - (52) 

^=1 ri7^?i,m£ „ — 1>/3 

From this we conclude by induction on /3 = M„ — 1, . . . , that 

x;:-^ = 0, if /3 < z/ < m^,„ - 1. (53) 
For /3 = M„ — 1, equation fl52|) reads 

n(C) L„ 

and thus (l53l) holds for f3 = Mn by the linear independence of the vectors 
Vq" (cf. Definition 12.51 (iv)). For the induction step we use (l52l) with (3 — 1 
instead of /3. Together with (153|) we find 



n(C) L„ 

^/3-l' 



0= n^^--^")"'" E ^""^ 
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which shows that flS51) holds for /3 — 1. Thus we have shown x = and this 
finishes the proof. □ 

The computational procedure is now a straightforward generalization of 
Section O First compute Bo,Bi e C^™'^' from (HE]). In addition to flM]) . 
assume 

mnk{W^V ■■■ A^-^W"V) = k. (54) 

Let us abbreviate 



/ V \ 



Compute the SVD 



%]iyr|] =Bo = K)So< 



H 



where Vq E C^'"'^ V^^^Vq = 4, Sq = diag((Ti, . . . , fXfc) G C'='^ and e C^'-'^ 
M^o^W^o = 4- From the rank conditions (l5Ti) .(l& 



(Ji > . . . (Tfc > = (Tfc+i = . . . = CTxi. 

The rank condition (ISTIl also implies 

i?(5o) = = i?(yo). 

Thus the matrix S = V^^V[K] G C^''^ is nonsingular and satisfies 

V[K] = VoS. 

With (go]), ([55]) we find 

W^Ij = S-'EoWo", 

and then from fl50|) 



(55) 



Finally, this leads to 

D := Vq^B^Wo^o^ = SAS-\ 
Therefore, the analog of Theorem 13.31 is 



(56) 
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Theorem 5.2. Suppose that T G H{Q, C"''") has no eigenvalues on the 
contour F in and pairwise distinct eigenvalues A^j-n. = 1, . . . ,n{r) inside 
r with partial multiplicities mi^„ > • • • > 'n^L„,n,'^ = 1, • • • ,"^(1^). Assume 
that the rank conditions (l5Ti) . (lMI) are satisfied with k given by (|T5|) . Then 
the matrix D G C'^''^ from (l56l) has Jordan normal form (l29i) with the same 
eigenvalues A„ and partial multiplicities „ = 1, . . . , L„, n = 1, . . . , n(r)). 
Suitable CSGEs for T can be obtained from corresponding CSGEs s^'^ for D 
via 

= ^J'^sf' ^<J< ^i,n - 1, 1 < ^ < ^n, 1 < < n(r), 
where Vq^^ is the upper m x k block in 



Vo 



(57) 



Remark 5.3. In a sense this generalization is similar to linearizing a poly- 
nomial eigenvalue problem by increasing the dimension. Note, however, that 
this only becomes necessary if there are too many eigenvalues inside the 
contour, or if rank defects occur that are not present in linear eigenvalue 
problems. 



The generalization of the algorithm from Section 13.31 is the following. 
Integral algorithm 2 

Step 1: Choose numbers I < m, K > 1 and a matrix V G C™'' at random. 
If more than m eigenvalues are expected inside F, let / = m, = 

Step 2: Compute 

and form i?o,Ar,-Bi,Ar as in fHHj) . 

Step 3: Compute the SVD Bq^n = VHW" , where 

V G C^'".^', W G C-^''^', V^V = W^W = Iku S = diag(ai, aa, . . . , am). 

Step 4: Perform a rank test for S, i.e. find < A; < Kl such that 

(Ji > . . . > o-fc > o-fe+i . . . am 

\i k = Kl then increase / or K and go to Step 1. 

Else let Vo = V{1 : Km, l:k),Wo = W{1 -.Kl,!: k) and 

So = diag(ai,o-2, . . . ,(1^). 

Step 5: Compute D = Vq^Bi^nWq^o^ G C^'^ 
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Step 6: Solve the eigenvalue problem for D 
DS = SA, S = {si. . .Sk),A = diag(Ai, . . . , A^). 

If ||T(Aj)fj|| is small and Xj G int(r) accept Vj = Vq^^Sj (with Vq^^ from (137|) ) 
as eigenvector and Xj as eigenvalue. 

5.2 Numerical Examples 

Example 5.4. We apply the integral algorithm 2 to the rank deficient ex- 
ample (H6l) . where K = 2,1 = 3 and the contour is the circle from fH5l) . Now 
the eigenvalues a = —0.2 and 6 = 1 are reproduced correctly (see Figure 
[6](left)), and three singular values survive as expected (Figure [6] (right)). 




Figure 6: Example 15.41 Eigenvalues from polyeig (open circles) and eigen- 
values from the integral algorithm 2 {K = 2, filled boxes) for a quadratic 
matrix polynomial with rank defect (left), singular values of integral algo- 
rithm 2 with / = 3 columns versus the number of quadrature nodes for 
the same example (right). 



Example 5.5. Consider the characteristic equation of a delay system x = 
Tox{t) + Tix{t - r) from [16, Sec.2.4.2],[llJ, given by 

T{z) = zI-To-Tie-^\ = ("2^ -g) ' = {~A -l) ' ^^^^ 

In case r = 1 there are more than two eigenvalues inside the circle ip{t) = 
zq + Re^\ ^ = —1, R = Q . We set Z = 2, = I2 and = 3 for the integral 
algorithm 2 and obtain with = 150 five eigenvalues inside the circle, (see 
Figure [Tl^left)), which coincide with the computed ones in [11]. Much smaller 
values than = 150 give sufficient accuracy, since there is a good separation 
of singular values and a fast decay of residuals, see Figure [7]( right). 
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Figure 7: Example [5]5l Eigenvalues of the characteric equation f l58|l inside a 
circle of radius 6 and with center —1, computed with the integral algorithm 
2 with K = 3,1 = 2. (left), residuals \\T{Xj)vj\\ for Ai ^ -0.6 + 2.71i, 
A2 ~ —2.27 + 5.07i versus the number N of quadrature nodes for the same 
example (right). 



References 

[1] T. Betcke, N. J. Higham, V. Mehrmann, C. Schroder, and F. Tisseur. 
NLEVP: A collection of nonlinear eigenvalue problems. Tech- 
nical Report 2008.40, MIMS, University of Manchester, Apr. 2008. 
www . mims . manchester . ac . uk/research/ numer ical-analysis/nlevp . html. 

[2] T. Betcke and H. Voss. A Jacobi-Davidson type projection method for 
nonlinear eigenvalue problems. Future Gener. Comput. Syst., 20:363- 
372, 2004. 

[3] P. J. Davis. On the numerical integration of periodic analytic functions. 
In On numerical approximation. Proceedings of a Symposium, Madison, 
April 21-23, 1958, Edited by R. E. Langer. Publication no. 1 of the 
Mathematics Research Center, U.S. Army, the University of Wisconsin, 
pages 45-59. The University of Wisconsin Press, Madison, 1959. 

[4] P. J. Davis and P. Rabinowitz. Methods of numerical integration. Dover 
Publications Inc., Mineola, NY, 2007. Corrected reprint of the second 
(1984) edition. 



29 



[5] I. C. Gohberg and E. I. Sigal. An operator generalization of the log- 
arithmic residue theorem and Rouche's theorem. Mat. Sb. (N.S.), 
84(126) :607-629, 1971. 

[6] R. E. Greene and S. G. Krantz. Function theory of one complex variable, 
volume 40 of Graduate Studies in Mathematics. American Mathematical 
Society, Providence, RI, third edition, 2006. 

[7] N. Hale, N. J. Higham, and L. Trefethen. Computing A"', log{A), and 
related matrix functions by contour integrals. SIAM J. Numer. Anal., 
46:2505-2523, 2008. 

[8] N. J. Higham. Functions of Matrices. SIAM, 2008. 

[9] M. V. Keldysh. On the characteristic values and characteristic functions 
of certain classes of non-self- adjoint equations. Doklady Akad. Nauk 
SSSR (N.S.), 77:11-14, 1951. 

[10] M. V. Keldysh. The completeness of eigenf unctions of certain classes 
of nonself adjoint linear operators. Uspehi Mat. Nauk, 26(4(160)):15-41, 
1971. 

[11] D. Kressner. A block Newton method for nonlinear eigenvalue problems. 
Numer. Math., 114:355-372, 2009. 

[12] A. S. Markus and E. I. Sigal. The multiplicity of the characteristic 
number of an analytic operator function. Mat. Issled., 5 (3 (17)): 129- 
147, 1970. 

[13] V. Mehrmann and H. Voss. Nonlinear eigenvalue problems: a challenge 
for modern eigenvalue methods. GAMM Mitteilungen, 27, 2004. 

[14] R. Mennicken and M. Moller. Root functions, eigenvectors, associated 
vectors and the inverse of a holomorphic operator function. Arch. Math. 
(Basel), 42(5):455-463, 1984. 

[15] R. Mennicken and M. Moller. Non-self-adjoint boundary eigenvalue prob- 
lems, volume 192 of North-Holland Mathematics Studies. North-Holland 
Publishing Co., Amsterdam, 2003. 

[16] W. Michiels and S.-l. Niculescu. Stability and Stabilization of Time- 
delay Systems, volume 12 of Advances in Design and Gontrol. Society 
for Industrial and Applied Mathematics (SIAM), 2007. 



30 



[17] S. Solovev. Preconditioned iterative methods for a class of nonlinear 
eigenvalue problems. Linear Algebra AppL, 415:210-229, 2006. 

[18] G. Stewart and J. G. Sun. Matrix Perturbation Theory. Academic Press 
Inc., Boston, MA, 1990. 

[19] V. P. Trofimov. The root subspaces of operators that depend analytically 
on a parameter. Mat. Issled., 3(vyp. 3 (9)):117-125, 1968. 

[20] H. Voss. A maxmin principle for nonlinear eigenvalue problems with 
application to a rational spectral problem in fluid-solid vibration. Appl. 
Math, 48:607-622, 2003. 

[21] H. Voss. An Arnoldi method for nonlinear eigenvalue problems. BIT, 
44:387-401, 2004. 

[22] H. Voss. A Jacobi-Davidson method for nonlinear and nonsymmetric 
eigenvalue problems. Comput. Struct., 85:1284-1292, 2007. 

[23] H. Voss and B. Werner. A minmax principle for nonlinear eigenvalue 
problems with applications to nonoverdamped systems. Math. Methods 
Appl. Sci, 4:415-424, 1982. 



31 



